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Abstract 

We propose a new numerical technique for following the evolution of a self-gravitating col- 
lisionless system in general relativity. Matter is modeled as a scalar field obeying the coupled 
Klein-Gordon and Einstein equations. A phase space distribution function, constructed using 
covariant coherent states, obeys the relativistic Vlasov equation provided the de Broglie wave- 
length for the field is very much smaller than the scales of interest. We illustrate the method 
by solving for the evolution of a system of particles in a static, plane-symmetric, background 
spacetime. 
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1 Introduction 



There now exists techniques to calculate numerically the dynamical evolution of matter in 
general relativity. For the most part research has focused on three very different types of 
systems; fluids, scalar fields, and collisionless matter. In this work we develop and exploit a 
connection between the latter two. 

A collisionless system with a large number of particles is treated generally as a continuous 
fluid in phase space. The fluid is described by a distribution function / which gives, for each 
region of the system, the density of particles and their distribution in velocity space. / obeys 
the Vlasov equation with an appropriate force law which, in general relativity, is given by the 
geodesic and Einstein equations. Shapiro and Teukolsky [1] have developed a computational 
method for handling collisionless relativistic systems that combines N-body techniques and 
numerical relativity. An alternative approach [2] is to follow the evolution of / directly in 
phase space by solving the coupled Einstein and Vlasov equations. Applications include violent 
relaxation and the collapse to a black hole of an unstable relativistic star cluster. 

In this work we show that a massive scalar field obeying the coupled Klein-Gordon and 
Einstein equations provides an alternative model for collisionless relativistic systems that can 
be readily adapted to the computer. The technique makes use of what is essentially the scalar 
field analogue of geometric optics in general relativity [3] . The spirit and methodology is similar 
to that of Widrow and Kaiser [4] who showed that a field obeying the coupled Schrodingcr and 
Poisson equations could be used to model collisionless, nonrelativistic matter. 

Scalar fields in general relativity have been studied for some time, albeit for entirely different 
reasons. For the past fifteen years the motivation for much of this work has been to understand 
the inflationary universe paradigm. For example, numerical simulations of inhomogeneous scalar 
field cosmologies have been studied in an attempt to understand the onset of inflation [5]. 
Recently, there has been great deal of interest in the gravitational collapse to a black hole of 
a self-gravitating massless scalar field. This is due primarily to the discovery by Choptuik [6] 
that such systems exhibit scaling behaviour and critical phenomena. 

Again our interest here is in simulations of collisionless matter. For our purposes previous 
investigations of scalar fields are important in that they exhibit the numerical techniques used to 
follow their dynamics. Our discussion does suggest that the results found by Choptuik might also 
apply to collisionless matter. Unfortunately the simulations that have shown scaling behaviour 
and critical phenomena involve massless scalar fields whereas our work requires a nonzero mass. 

We model a collisionless relativistic system by a field configuration cf)(x) where is a complex 
scalar field obeying the Klein-Gordon equation in curved space: 



(The metric has signature (+ — ) and we set c = 1.) m/Ti is a model parameter which 

must be large enough to guarantee that geometric optics (or rather the scalar field analogue 
of geometric optics) applies. In particular we require that throughout the system A dcB <C C 
and AdcB "C TZ where £ is the scale over which the density and velocity field of the system 
(i.e., distribution function) vary, 1Z is the typical radius of curvature for the spacetimc, and 
AdcB ~ \<f>/V<f>\ is the typical de Broglie wavelength for the field. Locally <j> can be regarded as 
the superposition of plane waves propagating along geodesies [3] . The amplitude of the field in 
a given region of the system tells us something about the density of particles while information 
about the velocity distribution of the particles is encoded in the phase structure of the field. 

The discussion above suggests the following prescription for constructing a distribution func- 
tion T = T(x,i>) from the field <fi: Fourier transform the field in a region of size r\ (AdeB <C rj and 
T) <C £, TV) centered on a particular point in phase space and then take the absolute square. This 
is what is done in [4] for a nonrelativistic Schrodinger field ip = ip(x, t). There T is constructed 
using coherent states: ,F(x, p) = |(e (x, p) |^)| 2 where |e (x, p)) is a localized state centered 
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on x with average momentum p. |e (x, p)) is taken to be a Gaussian (minimum uncertainty) 
wavepacket, 

/ , x 3/2 / -. v 3/4 

(e (x, P ) |x'} = ( -U -U /V-^'.p/n ; (2) 



y 27T?l / \ 7T77 2 / 

having width 77 in position space and hjr\ in momentum space. The normalization is chosen so 
that J d 3 x J d 3 pj r (x,p) — 1 provided J d 3 x\ip(x)\ 2 = 1. It is straightforward to show that T 
obeys approximately the nonrelativistic Vlasov equation provided AdcB 

The coherent-state representation for a scalar field in special relativity is again based on a 
set of state vectors \e(x, p)) localized in both position and momentum. We choose minimum 
uncertainty wavepackets [7]: 

/ 2 \ 3/4 

(e(x,p)\p') = (ij e-' iy/S - p '» J/s! (3) 

where p = (po, p) with po > |p| and a ■ b = ao&o — a • b. These wavepackets are centered about 
x and move with average momentum p. 

The extension to general relativity is most easily accomplished using Riemann normal coor- 
dinates (RNCs) and is similar to discussions found in the literature [8] of the Wigner func- 
tion in curved space. In a RNC system curved space closely resembles flat space in the 
neighborhood of a particular point taken to be the origin. The metric can then be writ- 
ten as a Taylor series about the flat space metric rj^ in increasing orders of the Riemann 
tensor R^ U p a '- g^v = f)\*u + \Rp.pvay p y a + • • • • In RNCs the phase factor x ■ p — 
rj^x^p 11 + 'post — geometric optics corrections [3]. 4> in the curved space coherent-state repre- 
sentation is then 

(e (x, p) \<j>) = j dV p e~ ix - p'/k-vp'v 2 /^ ^ (4) 

where 

_ dp dp 1 dp 2 dp 3 , 2 2 s 

dV P = 77~4 2lT 5 (P ~ 171 ) e (Po) ( 5 ) 

(2tt) 

is the momentum space volume element on the mass shell, Q(x) — 1 for x > and zero otherwise, 
and 4>(p) is the Fourier transform of <^(a;): 



dV p e- ip ^ h $(p) . (6) 

Note that (e (x, p) \<j>) reduces to (e (x, p) \ip) for |p| <C m. In practice, one integrates over po 
making use of the identity J dV p = d 3 p/ ((27r) 3 2p°). For example 

F{x,p) = f dV p 'dV;e* x i p '- p '') /h ~ p < p ' +p '') v2/h2 4>{p')4>*(p'') (7) 
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d 4 p' d 3 p' 
(2tt) 3 2p'° (2tt) 3 2p"° 



e *x{p'-p'')/fr-p{p'+p'')ri 2 /K 2 ( j ) (j ) >)j ) *( p >>y ( 8 ) 



In the second of these expression the integrand is evaluated "on mass-shell", i.e., with p = 
-^/p' • p' + m 2 and p"0 — y/ p" ■ p" + m 2 . 

Consider the function a(p) = S (p 2 — m 2 ) 9 (po) 4>{p)- From Eq. ([!]) and the definition of the 
RNCs we have 

1 " r " 2 \ R V«iJ^rP^ ~ Wp-Stp*) "00 - o(p) ff (p) = (9) 



h 2 3 p ^ dppdp(j 3 "dp, 



p 
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from which follows immediately the identity 



= [ d 4 p'd 4 p"e tx i p '- p ") /h - p < p ' +p ") v2/R2 (0(p')a(p'))a*(p") . (10) 



A straightforward but tedious calculation leads to a complex equation for T which can be split 
into real and imaginary parts: 



2 2 

vnr - p 



T = (11) 



The ". . ." in these equations refer to terms higher order in the rj/L, XdeB/v and rj/TZ (i.e. post- 
geometric optics corrections) . The first of these equations tells us that the distribution function 
is peaked on the p^p^ = m 2 (mass shell) surface. The second equation is, to leading order, the 
relativistic Vlasov equation in a RNC system. 

For self-gravitating collisionlcss matter, g^ depends on / through the Einstein equations, 
Rpv — \g^vR = 87rGT M „, where the stress-energy tensor T M1/ is given by 

[ d 4 p 

T^v = I —=p^p v f{x,p) . (13) 



We can calculate T pt/ for our model by substituting T for / in Eq. (13). This is, however, 
computationally expensive (one must update g at each timestep) and not really necessary. 
Consider instead the expression 



T„ v = d^d^* (14) 
= J dV p ,dV p/ >e i ^'-^ x /%p , J{p , )4,*{p'') . (15) 

Eq. (|l3|) is what one would obtain if this expression were smoothed over a length scale ~ r\ and 
since the quantity of interest, g ~ V -2 T, is itself smoother than T we can use either expression. 
Eq.(12) is just the usual expression for the stress energy tensor of a scalar field, save for a term 
proportional to d^d^cjf — m 2 <jxj)* . But this term is negligible in the geometric optics limits [9]. 

Evidently a self-gravitating scalar field behaves like collisionless matter and can be repre- 
sented in phase space using the coherent-state representation described above provided Idb <C 
r\ -C C, 1Z. Of course to follow the field properly we must choose a grid spacing that is somewhat 
less than Ade b ■ In addition the timestep must be short enough to follow the rapid evolution of 
the phase factor for <p. The accuracy of a simulation is limited by constraints on memory and 
CPU time. The situation is similar to what is encountered in particle methods where a finite 
number of particles is used to provide a statistical description of the distribution function. One 
can improve the accuracy of a given simulation by increasing the number of particles but again 
this comes at the cost of memory and CPU time. 

As an illustration of our method we consider "particles" in a static and plane symmetric 
background spacetime. For simplicity, we choose the line element ds 2 — A 2 dt 2 — A~ 2 dz 2 with 
A 2 = exp (lo 2 z 2 ) so that nonrelativistic particles near the origin will execute simple harmonic 
motion. It is convenient to write our equations in terms of the dimensionless quantities £ = loz, 
t = tot, k = p/m and C = m/huj. Eq. (Q), for example, becomes 
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AdcB should be made as small as possible (this is done by increasing L) but no smaller than 
several grid spacings. With this in mind, we choose C ~ 0.1 AT and -q ~ £ 1/<2 where N is the 
number of gridpoints used in the simulation. There are roughly TV 1 ' 2 resolution elements in 
both position and momentum space which is roughly what one would expect for an N-body 
simulation with TV particles. 

Consider first an initially cold (zero velocity dispersion) distribution of particles. The field 
configuration at r = is taken to be <f)(() = exp (— C 2 /Co) an d d<f>/dt — —i&fi with Co = 0.5. 
The corresponding T is shown in the top left panel of Fig 1. The phase space configuration at 
t = An (i.e., time when particles near the center of the distribution have made two complete 
orbits) is shown in the bottom left panel. For comparison, the results of an N-body calculation 
are shown in the right-hand panels. As expected, particles near the center of the distribution 
execute simple harmonic motion while particles initially at z > 1 follow anharmonic, relativistic 
orbits in phase space. 

We next consider a system of "hot" particles. In particular, we assume the system is initially 
described by a truncated isothermal distribution function: /(e) = 0(e) exp (e/T") where e = 
Po/E max and T" = T/E max . For our simulation, we take -E max = 2m and T = m. po is 
conserved for each particle in the system and therefore /(e) is constant. 

The initial wavefunction for the Klein-Gordon simulation is taken to be 

1 N 

^ = JjT.ieixuPiMx) (17) 
i=i 

where the pairs (xi,Pi) are chosen at random from the distribution function. This is reminiscent 
of N-body simulations where a statistical representation of the initial distribution function is 
constructed by specifying the positions and velocities of N super-particles. Much of the art in 
N-body simulations is in calculating the forces and various methods (e.g., particle-particle with 
smoothing, particle-mesh) are used to deal with such problems as artificial two-body relaxation. 
Here, the discrete particles are replaced by wavepackets and so smoothing is built in. The pros 
and cons of our method, as compared with N-body techniques, will be explored in a future 
publication [10]. 

Fig. 2 shows the distribution function T for r = and r = 87r; Fig. 3 gives !F(po) (calculated 
from !F(x 7 p) using a simple binning procedure) again for r = and r = 87r. While individual 
features of T change (this also occurs in particle realizations of /) the overall structure of the 
distribution function remains constant. 

This work illustrates that simulations of collisionlcss matter in general relativity can be done 
by following the dynamics of a massive scalar field and in so doing, bridges the gap between 
two very separate branches of numerical relativity. Applications to self-gravitating systems with 
more complicated geometries are straightforward only because of the large volume of work on 
relativistic scalar fields that already exists. 
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Figure Captions 

FIG.l. Evolution in phase space of an initially cold distribution of particles. The system is 
set in a fixed background spacetime as described in the text. 4096 grid points are used for the 
Klein-Gordon simulation. Left-hand panels give the Klein-Gordon distribution function T while 
right-hand panels give the corresponding N-body distribution function /. 

FIG. 2. Evolution of a "hot" system in phase space. The initial distribution function is that 
of a truncated isothermal sphere and is shown in the left panel. The distribution function for 
t = Stt is shown in the right panel. 2048 grip points are used for the simulation. 

FIG. 3. Distribution function T(E). T{E) is calculated from the T{x,p) of Fig. 2 using a simple 
binning procedure. Solid curve is the initial F(E); dashed curve is the T(E) at r = 87r; dotted 
curve is the F{E) used to construct the initial wavefunction. 
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